home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / general / fractal / kaos.lha / complib / nl_interpol.c < prev    next >
Encoding:
C/C++ Source or Header  |  1989-11-18  |  1.6 KB  |  63 lines
  1. /*
  2. ### non-linear interpolation by iteration of Poincare section ###
  3.  
  4. Note:     vx, vx2:oin primary coordinates
  5.     vseci1,vseci2: in section coordinates
  6. */
  7.  
  8. nl_interpol(tso,vx2,vseci1,vseci2,vx,time,tsi,dim)
  9. double vx2[],vseci1[],vseci2[],*tso,tsi,time;
  10. int dim;
  11. {
  12.     int i,it,it_max;
  13.     double ts,ts1,ts2,nl_interpol_f(),eps,err,fabs(),*vsec1,*vsec2,*vsec,*dvector();
  14.     void free_dvector();
  15.     extern int stop,nok,nbad,full_dim,polar_section,section_index;
  16.     extern int (*f_p)();
  17.     extern double section_constant,*t_vf,*t_v;
  18.  
  19.     vsec = dvector(0,full_dim-1);
  20.     vsec1 = dvector(0,full_dim-1);
  21.     vsec2 = dvector(0,full_dim-1);
  22.     eps = 1.e-14;
  23.     it_max = 10;
  24.  
  25.     ts1 = 0;
  26.     ts2 = tsi;
  27.     for(i=0;i<full_dim;i++) vsec1[i]=vseci1[i];
  28.     for(i=0;i<full_dim;i++) vsec2[i]=vseci2[i];
  29.  
  30.     ts = ts1 + (ts2-ts1) * (section_constant-vsec1[section_index])/(vsec2[section_index]-vsec1[section_index]);
  31.  
  32.     for(it = 0;it <it_max;it++){
  33.         int_onestep(t_v,t_vf,vx,&time,ts,dim);
  34.         to_window_variables(vsec,t_v,polar_section);
  35.         if((vsec[section_index]-section_constant)*(vsec1[section_index]-section_constant)) {
  36.             ts1 = ts;
  37.             for(i=0;i<full_dim;i++) vsec1[i]=vsec[i];
  38.         }
  39.         else {
  40.             ts2 = ts;
  41.             for(i=0;i<full_dim;i++) vsec2[i]=vsec[i];
  42.         }
  43.         if(stop)
  44.             goto done;
  45.         if(fabs(ts-ts1)<eps)
  46.             goto done;
  47.         ts = ts1 + (ts2-ts1) * (section_constant-vsec1[section_index])/(vsec2[section_index]-vsec1[section_index]);
  48.     }
  49.  
  50. done:
  51.     *tso = ts;
  52.     err = fabs(ts-ts1);
  53.     from_window_variables(vx2,vsec,polar_section);
  54.     if(it>=10)
  55.         nbad++;
  56.     else
  57.         nok++;
  58.  
  59.     free_dvector(vsec,0,full_dim-1);
  60.     free_dvector(vsec1,0,full_dim-1);
  61.     free_dvector(vsec2,0,full_dim-1);
  62. }
  63.